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Abstract. A first- and second-order relation between cos- 
mic density and peculiar-velocity fields is presented. The 
calculation is purely Lagrangian and it is derived using the 
second-order solutions of the Lagrange-Newton system ob- 
tained by Buchert & Ehlers. The procedure is applied to 
two particular solutions given generic initial conditions. 
In this approach, the continuity equation yields a rela- 
tion between the over-density and peculiar-velocity fields 
that automatically satisfies Euler's equation because the 
orbits are derived from the Lagrange-Newton system. This 
scheme generalizes some results obtained by Nusser et al. 
(1991) in the context of the Zel'dovich approximation. As 
opposed to several other reconstruction schemes, in this 
approach it is not necessary to truncate the expansion of 
the Jacobian given by the continuity equation in order to 
calculate a first- or second-order expression for the den- 
sity field. In these previous schemes, the density contrast 
given by (a) the continuity equation and (b) Euler's equa- 
tion are mutually incompatible. This inconsistency arises 
as a consequence of an improper handling of Lagrangian 
and Eulerian coordinates in the analysis. Here, we take 
into account the fact that an exact calculation of the den- 
sity is feasible in the Lagrangian picture and therefore an 
accurate and consistent description is obtained. 
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1. Introduction 



An important aim of cosmology is to determine a rela- 
tionship between the cosmic density contrast d and the 
peculiar-velocity field u under the assumption that these 
fields have evolved under the action of gravity. In the Eu- 
lerian theory of gravitational instability, retaining only the 



linear growing-mode solution, this relationship is simply 
given b}|] 

S{t,q)^^^V,-u{t,q), (1) 
aD 

where D{t) is the usual linear growth factor, such that 
D{t) = a{t) = {t/to)'^/^ in an Einstein-de Sitter universe. 
Beyond linear theory, in the weakly non-linear regime, 
has been used to good effect as an estimate of the den- 
sity contrast field from observed peculiar-velocities (e.g. 
Bertschinger & Dckel 1989; Dekel et al. 1990; Bertschinger 
et al. 1990). On the other hand, the POTENT method and 
other quasi-linear reconstruction procedures involve the 
Zel'dovich approximation (ZA) (Nusser et al. 1991; Nusser 
& Dekel 1992). The ZA does not however, strictly speak- 
ing, provide an exact algebraic relation between the fields 
u{t, q) and S{t, q). It can however yield (as shown in Nusser 
et al. 1991) a relation between the fields that is self-consistent 
within the order of the approximation, i.e. within 0{D^). 
In the ZA, from the Lagrangian integral of the continuity 
equation we have 

-1 



l + 5{t,X) 



= det 



D dui{t,X) 



(2) 



aD -^.j 

and the fields are mapped back to Eulerian space via the 
ZA map for the orbits of the particles. Initial conditions 
are specified at an arbitrary early time tg > 0. Also, it is 
customary to make the assumption of smoothness at early 
times, (5(to, X) « 0, by virtue of the amplitude of the fluc- 
tuations in the microwave background. From (H) we obtain 
the linear relation (|^) by truncating the expansion of the 
determinant to the lowest order in the perturbation 
field u and replacing Lagrangian by Eulerian coordinates 
{X^q). 

These schemes have been generalized further (Gramann 
1993a, 1993b; for their Eulerian analogue, see Bernardeau 
1992) by starting out with a parametrization of the parti- 
cle orbits q — F{t, X) of the type: 

F{t,X) = X + Dy x'i>'^^HX) + D'^V x<^'-^\X), (3) 



^ Throughout this article, q and X denote (comoving 
lerian and Lagrangian coordinates respectively. 
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In the case of a flat Universe, the leading second-order 
term in perturbation theory (Bouchet et al. 1992) is in- 
deed oc D^{t). Therefore, the coordinate map (g) is in 
principle well-motivated (for a numerical comparison of 
the different models, see e.g. Dekel 1994). However, as 
shown in Gramann (1993b), orbits of the type (^) do not 
yield a self-consistent relation for the cosmic density and 
velocity fields; the density fields obtained via (a) the con- 
tinuity equation and (b) the Eulcr equation differ within 
the order of the approximation. The reason for this is that, 
as we shall see, (||) is only an approximate second-order 
solution of the Lagrange-Newton system. Thus, we wish 
to call the attention of the reader to the fact that the re- 
construction models derived from and, more generally, 
from an ansatz including higher orders in D(t), 

N 

F(t,X) = X + ^L>" Vx$'"HX) (4) 

(for arbitrary N), do not follow a rigorous line of analysis 
(in the case of an arbitrary il) for the following reasons: 

I In these models, the parametrization of the orbits is 
given ad hoc and it is not derived as a solution of the 
Lagrangian evolution equations for the flow field. In 
particular, the absence of lower-order growing modes, 
that are present in the general perturbation theory so- 
lution, implies that the derivation of d through the 
continuity equation neglects couplings of these modes 
with leading-order terms. 

II The density contrast Sc obtained from the continuity 
equation in perturbation theory and Se satisfying Eu- 
ler's equation are mutually incompatible. 

III There is no reason to justify that D{t) is a good per- 
turbative parameter (in terms of the convergence of 
the solutions) and that a meaningful generalization of 
the family of solutions where the ZA belongs can be 
realized as a polynomial in D(t). 

It is easy to see that (II) follows as a consequence of (I). 
Regarding (I) and (III), it has been shown that the cor- 
rect coordinate map between Eulerian and Lagrangian co- 
ordinates in Lagrangian perturbation theory is obtained 
by calculating the perturbativc solution of the Lagrange- 
Newton system for the trajectories q = F{t,X). These 
evolution equations are obtained by transforming the Euler- 
Newton system to Lagrangian coordinates and eliminating 
all Eulerian fields by using exact Lagrangian integrals for 
the Eulerian acceleration and density (Buchert & Gotz 
1987; Buchert 1989). Within this scheme, the ZA is re- 
covered as a subclass of solutions (Buchert 1992), and 
hence, one can obtain a self-consistent solution for the 
over-density field in terms of the velocities as in Nusser et 
al. (1991). Furthermore, the second- and third-order so- 
lutions obtained in the Lagrangian framework (Buchert 
& Ehlers 1993, Buchert 1994; and in a shghtly different 
approach Bouchet et al. 1992, 1995; Catelan 1995) are 
unique and well-defined at all times from a fiducial time 



to, where initial data are specified,^ until shell-crossing, 
provided we impose periodic boundary conditions and fix 
some global gauge conditions (Ehlers & Buchert 1997) as 
opposed to the orbits described by polynomials in D{t), 
which are expected to have poor convergence (for a dis- 
cussion in the context of the variational approach to the 
solutions see Susperregi & Binney 1994). Although the 
leading terms of the longitudinal parts of the perturba- 
tion solutions seem to confirm this polynomial approach 
(e.g. Bouchet et al. 1992, for the second-order solutions), 
this is not a mathematically consistent motivation for us- 
ing D{t) as an expansion parameter by going to higher 
orders in perturbation theory. On the contrary, the struc- 
ture of the Lagrangian perturbation scheme is such that, 
starting at third order, the presence of interaction terms 
among perturbations of different orders are an obstacle to 
convergence, and in fact for ti > 3 all terms are interaction 
terms (see Buchert 1994). Therefore, for n > 3 an ansatz 
of the form (^ cannot reproduce gravitational dynamics, 
as stated in (III). This is an especially critical point to 
bear in mind in the context of reconstruction models in 
which higher-order couplings of flow fleld gradients show 
up in the density-velocity relation. 

In this paper, we approach the reconstruction problem 
from a purely Lagrangian perspective using the solutions 
of the Lagrange-Newton system obtained by Buchert & 
Ehlers (1993). The goal of the paper is to obtain a self- 
consistent solution for the over-density field that is valid 
to second order, thus extending the results of Nusser et 
al. 1991 in the context of the ZA. Our results are purely 
Lagrangian, but in obtaining them we rely on the fact that 
there is a one-to-one mapping between Lagrangian and 
Eulerian coordinates and hence, the reconstruction for- 
mulae presented are only valid up to orbit-crossing time. 
The basic formalism is laid out in Sect. 2, and in Sect. 3 
we discuss its application to two different classes of ini- 
tial conditions. In Sect. 4 we discuss the self-consistency 
of solutions in perturbation theory models. Sect. 5 com- 
pares the model presented with previous reconstruction 
procedures and Sect. 6 sums up. 

2. First- and second-order Lagrangian field 
reconstruction 

In what follows, we shall mainly discuss the physically 
interesting case where the peculiar-velocity field is par- 
allel to the peculiar-acceleration, up to decaying modes, 
a property which results from the Eulerian linear theory 
for irrotational fields, and is also valid for the general ir- 
rotational Lagrangian linear solution (the growing-mode 
solution supports the tendency to becoming parallel, see 
Buchert (1992)). We will further recall that both fields are 

We make use of the term 'initial' to denote conditions given 
at to, an arbitrary but sufficiently early epoch, not necessarily 
t = 0. 
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practically parallel up to orbit-crossing time, also in the 
second-order schemes we will consider. The results pre- 
sented however are not restricted to solutions where this 
property holds but, as we will see below, apply to a wider 
class of solutions. 

In what follows we shall study subclasses of the second- 
order Lagrangian perturbation solutions given by Buchert 
& Ehlers (1993) which are only restricted by the require- 
ment that there exists an arbitrary functional relationship 
among peculiar-velocity and -acceleration at the initial 
time. The subclasses are singled out by special assump- 
tions for this functional relationship. 

The second-order parametrization of the orbits of the 
particles in terms of their Lagrangian coordinates and 
time for the subclasses of solutions which we consider later 
can be written as follows: 

q - F{t, X) = X+h (i) Vjf 5(1) {X)+h (i) Vx5(2) (X), (5) 

where fi{t) and f2{t) are dimensionless functions of time, 
calculated from solving the Lagrange-Newton system for 
deviations from the homogeneous solution to second order. 
Moutarde et al. (1991) investigated the orbits (^ (rather 
than (Q)), for leading orders of /i(i) and /2(t) and for a 
model of a particular initial potential. In the next section, 
we will examine two particular solutions for /i {t) and /2 {t) 
for different, but generic, initial conditions. 

The expansion factor a{t) is normalized to unity at to, 

a{t) — (j^) 1 (and equal to the linear growth factor D{t) 
in an Einstcin-de Sitter universe considered here). The 
perturbation potentials S*^^^ and S*'^' are determined by 
the initial data such that, in general, they have to satisfy 
the (non-local) boundary conditions: 

V|5«=/(5)io, (6) 
5(2) ^ 2 7/(^(1)), (7) 

where S denotes the initial peculiar-velocity potential 
u{tQ^X) = \JxS{X)^ and (following the summation con- 
vention) 



0(1) cW\ 



(8) 
(9) 



are the first and second principal invariants of the peculiar- 
velocity tensor gradient. For periodic boundary conditions 
and suitable gauge conditions, however, which we shall 
assume hereafter, we can simplify (||)(0) (see Ehlers & 
Buchert (1997) for a rigourous proof): 

VxS^'^ = VStQ , (10) 

V^S-^^) = 2//(S')tg . (11) 

Therefore, (q) gives us the trajectory of any given particle 
by solving (10|)(pT|) for an initial potential S. At the time 
to, both fi and /2 vanish, so Lagrangian and Eulerian 
coordinates coincide. 



Thus, the expression for the density contrast at any given 
time and for any trajectory field reads: 



[l + S{t,X)] = [l + 5{to,X)] J-i 
where the determinant J{t, X) is 



J = det 



dXi 



(13) 



(14) 



In order to calculate (5 as a function of Eulerian coordi- 
nates, we need to invert the coordinate map (|^) so that 
Lagrangian coordinates are expressed in terms of the Eu- 
lerian coordinates and thereby derive the following: 



[l + S{t,H{t,q))] = [1 + S] J' 



(15) 



where S{X) is the initial density contrast and the deter- 
minant J^{t,q) = J-^{t,X) is 



J 



^ - det - 



(16) 



q = F{to,X) = X. 



(12) 



Eq. (|T^) gives the Eulerian values of the field 5 at the 
Eulerian positions corresponding to a given trajectory la- 
belled hy X — H(t, q). In principle, we can interpret this 
as a field 6{t, q) satisfying the Eulerian equations of mo- 
tion. Nonetheless, it is important to bear in mind that 
the functional dependence oi 5 on H is only well-defined 
while the map between both coordinate spaces remains 
one-to-one. The peculiar-velocity is given by 

u{t,H) = act -^F , (17) 
da 

and, similarly, the peculiar-acceleration, 
d d^ 

w{t,H) ^ {aii + 2a^)—F + acP'—rF . (18) 
da da^ 

The Jacobian of the transformation H is given ex- 
plicitly in Appendix B. This general expression involves a 
number of cross-terms of u and w components. However, 
we can ignore these by assuming the peculiar- velocity field 
to be parallel to the peculiar-acceleration, such that 

w{t,H)^u{t,H)r^. (19) 

This assumption can be justified on physical grounds: after 
some time, the peculiar-velocity field tends to be parallel 
to the peculiar-acceleration, for the second-order irrota- 
tional solution (^), due to the existence of decaying and 
growing solutions in the weakly non-linear regime. We will 
address this issue in detail in Sect. 5, and show that a 
weaker condition on the gradients of the fields leads to 
the same result. As in the first-order model (where (19) 
holds exactly) the growing part supports parallelism, as 
it is shown in Buchert & Ehlers (1993) for a wide class of 
irrotational second-order solutions. Indeed, the difference 
between ut~^ and w for (||) is given by decaying modes 
only, as we shall see in the next section. Therefore, we are 
entitled to express w by some function proportional to u, 
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SO that the inversion of the map (|^) satisfies the implicit 
equation: 



H{t,q)^q~h{t) u{t,H{t,q)), 



(20) 



where h{t) is given explicitly in Appendix A for two dif- 
ferent models. 

In the first-order model, by virtue of (p^, the same 
relation ( |2C| ) holds with the appropriate h{t) (Appendix 
A.l). With condition ( p^ implicit in (pO[), the determinant 

reduces to: 

= l-h I^{u,) + Il^iu,) - e III^{ui). (21) 

Here, (divergence), 11^ (dispersion of diagonal compo- 
nents) and III^ (determinant) denote the three principal 
scalar invariants of the peculiar-velocity gradient with re- 
spect to comoving Eulerian coordinates as given in Ap- 
pendix B. Therefore, substituting ( [2l| ) into (|l^), we derive 
the relation for the density and velocity fields: 
l + 5{t,H{t,q)) _ 



1 + 5 

1 - h I^{ui) + Il^(ui) - in^{u,) 



(22) 



where Ui = Ui{t, H{t,q)). This result is unique and valid 
up to orbit-crossing time, i.e., while the coordinate inver- 
sion X = H{t, q) is well-defined. In principle, it enables us 
to calculate the cosmic density field from observations of 
the peculiar-velocity field in the weakly non-linear regime 
at a fixed time, provided we set 5{to) = 0. It is interest- 
ing to emphasize that, due to the parallelism condition 
(|l9|), ( p2|) applies at both first- and second-order (the ap- 
propriate h{t) for each case is given in Appendix A). The 
simplicity of ( p2| ) is remarkable. It shows that the first- 
and second-order models contain not only quadratic terms 
in Uij (as pointed out in the literature, e.g. Gramann 
1993a, b) but also cubic terms, arranged in the form of 
the principal invariants of the deformation tensor: I^{ui), 
II^{ui) and III^{u^). 

A more general result, albeit a complicated one, can 
also be derived for the case where we drop the condition 
(p^). The corresponding solution, taking into account all 
cross-terms, is obtained using the full solution for given 
in Appendix B. Nonetheless, considering w is not directly 
measurable from observations, the application of the gen- 
eral solution without the assumption of parallelism would 
be impractical. In these circumstances ( p2| ) provides a use- 
ful approximation. 



3.1. Model I: u{X) = w{X)to 

In this case, the velocity and acceleration fields are set to 
be identically parallel at the initial time. This relation is 
satisfied in Zel'dovich's ansatz, and it is a natural condi- 
tion to set in a more general scheme, given the remarkable 
accuracy of the "Zel'dovich approximation" . The initial 
density perturbation is thus expected to be proportional 
to the divergence of the peculiar- velocity. The solution ^ 
then is specified by the functions 



h{t) = -[a{t)-\], 



hit) , - . 

I 4' 14 ' ' 5 ' ' 2 35 

Therefore, the peculiar-velocity reads 



-3/2 



4- -( — a 
2^ 7 



3/2 



,1/2 



6 



5 35 
and the peculiar-acceleration, 
w{t,H) 



(2) 



(23) 
(24) 

(25) 

(26) 



Thus the departure from the parallelism condition (|9|) is 
3 



u t 



7 



(l-a-^/2)t-2 VxS^^\ 



(27) 



Hence, this difference is a negligible effect due to the pres- 
ence of the decaying mode. At i = toi the RHS of (27) 
is identically zero and, at later times, the decaying mode 
vanishes sufficiently fast and the overall coefficient there- 
fore tends to a constant. This constant second-order term 
could be absorbed in a coordinate transformation. In this 
case, h[t) in the reconstruction formula ( |2^ ) is given by 
(see Appendix A.2 for details) 

3 



■to (1 



i-V2)-i(ai/2 



7-1/2 



+ a-4). (28) 



This is a monotonic function of t, and to the leading order 
it reads 



(29) 



which is consistent with the results presented in Moutarde 
et ah (1991). 



3. Application to specific models 

In this section, we will examine two particular solutions 
of the second-order orbits (|^) that are of interest in the 
reconstruction procedure described in Sect. 2. These are 
characterized by the nature of the initial conditions at to ■ 
At this time, only one field, S, is specified. These solutions 
are derived in Buchert & Ehlers (1993). 



3.2. Model II: w{X) = Q 

Although Model I provides an appropriate initial value 
setting from the physical point of view, the initial density 
contrast in this model is strictly non-vanishing. However, 
for reconstruction models we have to neglect this, since we 
do not want to specify initial information. Precisely speak- 
ing, this implies an error which can be removed by using 
Model II: Given a vanishing initial peculiar-acceleration 
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field, the initial density contrast S is zero and its growth 
is induced by the initial velocity perturbations only. In 
this case, fi{t) and /2(i) in the class of solutions for the 
orbits (|^) are given by 



/iW = -[a(t)-a-3/2(t)], 



/2W = 



9 3 
25 ^^li'^ 



29 
40 



43 
— ( 
70 



-3/2 



it) 



(30) 



(31) 



Therefore, from (|l7j) the peculiar-velocity is 



u{t,H 
29 



1/2 



60 



2 

43 
70' 



o 1 



9 

25 ^ 



^3/2 



7 



(32) 



and the peculiar-acceleration. 



w{t,H) ^ - (a^' 
5 

3 , 10 29 1 
+ 25(-y+20" 



,-7/2w-2 



43 _7/9 5 _r;, 

35 4 ^ 



i-'Vx5(').(33) 



Like in the previous example we can calculate the depar- 
ture from parallelism which is given in terms of the decay- 
ing modes 

3 ,4 1 _5/2 43 _7/2 

H ( a a ^'^ - 

25^7 2 14 



such that, it is once more manifest that at late times par- 
allelism is supported by the growing modes. At first sight 
it might seem that the presence of the first-order poten- 
tial in (^^ is a shortcoming for model II in comparison 
to model I, where departure from parallelism is a second- 
order effect, given by (^^. However, although the initial 
fields are non-parallel as indicated by (|3^), the orders of 
magnitude involved in the RHS are very small from the 
outset and these decaying modes soon decrease further to 
negligible values. 

3. 3. Parallelism of u and w 

The assumption of parallelism ( p^ ) is justified by the pres- 
ence of decaying modes in the difference tw — u m both 
models I and II. On the other hand, it is apparent that the 
exact expression for requires knowledge of the deriva- 
tives of the peculiar-acceleration field. Using the paral- 
lelism condition to simplify the solution implies neglecting 
the contribution of the Eulerian derivatives of the vector 
field tw — u. In this section we produce an analytical argu- 
ment to justify this approximation in the weakly nonlinear 
regime. 



A measure of parallelism of u and w is given by the 
quantity |u x (an alternative measure has been pro- 
posed recently by Bagla & Padmanabhan (1996)). We ex- 
pect this quantity to be small in the regime where u and 
w are close to being parallel. As is given explicitly in Ap- 
pendix A, we have 

uit,H) ^a{t)\7xS''^^ + f3it)VxS^^\ (35) 

vmxS^^\ (36) 



and therefore, 

uxw=(aC-Pri) Vx5'(^) x VxS'^^^ 



(37) 



For the particular case of Model I, for instance, this leads 
to 

I uxw 1= ^t^^ I l-a-^/2 I |Vx^''^xVx5(2^|.(38) 

This quantity vanishes for all models where VxS"*^^^ and 
Vx'S''-^-' are parallel. However this is not so in general, 
and the integrals of (||)(@) can yield to non-parallel (and 
non-neglig ible) gradients VxS'^^^ and \7xS'^^\ Further- 
more, the presence of the coefficient a^^'^{t) in ( |38| ) in- 
dicates that an initial departure from parallelism of the 
fields u and w will grow in time. At this point it is impor- 
tant to stress that our simplification to derive (|2^) does 
not strictly speaking rely on the parallelism of the fields u 
and w but on that of the spatial gradients Uij and Wij . 
The differences of these fields are intimately connected, 
but the time dependence involved is different as is clear 
by direct differentiation. Parallelism of the gradients is a 
weaker condition to impose than that of the fields them- 
selves. It is easy to show that 



H,3 



4, +/l Vfc,^«+/2 Vfe,^(2) 



(39) 



and 



4-, +/i Vfe,5(i)+/2 Vfe,5(2 



Therefore, we have that for Model I: 



Wi. 



7(1 



(40) 



(41) 



and for Model II: 



r 



5 



kj 



25 V7 



- ia-5/2_ 15^-7/2 
2 14 2 



where 



V,fc5(2) rfc,,(42) 



(43) 
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In the case of Model I, parallelism is exact in the ini- 
tial conditions, therefore there is no dependence of the 
derivatives on the potential S'^^K In all other models, the 
difference Uij — Wij comes as a linear combination of 
the double derivatives V^jS''-"-', n — 1,2, up to the co- 
efficient Tij that plays a role when we approach orbit- 
crossing. For realistic fields, i.e., those that evolve from 
negligible amplitudes at t = Iq, the gradients of the per- 
turbative potentials S''^\ 5'^^'' are sufficiently smooth and 
slowly varying, and therefore they are well under the or- 
der of the approximation in the weakly nonlinear regime. 
We stress, however, as indicated above that in the strongly 
nonlinear regime (or at the epoch of shell-crossing, respec- 
tively) , this argument breaks down in the same manner as 
perturbative solutions at any order do. 



4. Self-consistency of perturbative solutions 



the orbits we have used are obtained by solving the Lagrange- 
Newton system perturbatively, Eulerian quantities con- 
structed from these will automatically satisfy all equa- 
tions to second-order. However, the density contrast con- 
structed from the velocities through ( p^ ) is associated to 
a Newtonian potential causing an acceleration that will in 
general differ from the acceleration given by the second- 
order orbits. This difference vanishes only for exact solu- 
tions. 

Let us denote by a (1) superscript the quantities de- 
rived from the particle orbits. Hence, 



- hit) i^iu\'>) + h^it) 11^ [uY') - in 

and therefore, 



In our reconstruction procedure we have proposed a sim- 
plification of the exact solution ( [l5| ) under the assump- 
tion of parallelism of the fields u and w. Thus we are able 
to obtain an expression for the density in terms of the 
derivatives of the velocity only. We have seen that this 
assumption is sufficiently accurate, the difference between 
both fields is a negligible effect. However, it is not nec- 
essary to adopt parallelism, and in this section we will 
outline a possible scheme to construct iteratively a "self- 
consistent" solution for w from the velocity field, so that 
the reconstruction method with the Jacobian (Bl) can be 
employed. 

At this point, we would like to call the attention of 
the reader to the problem of self-consistency of approx- 
imations as pointed out by Doroshkevich et al. (1973) 
(compare the discussion by Buchert (1989)). We enun- 
ciate the problem by stating that perturbative solutions 
of the equations of motion are mutually consistent only 
within the truncation order. This can be illustrated with 
the "Zel'dovich approximation" . The density contrast as 
given by the Poisson equation is 

3 



(44) 



and the density contrast as given by mass conservation 



J{t,H)- 

Introducing the dimensionless error 

_ _ ^(1) 



(45) 



(46) 



this difference for the "Zel'dovich approximation" gives 

e^{a~l)^ II{S) + {a-l)^ III{S), (47) 

which is within the order of the approximation. In the 
second-order model we expect an analogous situation. As 



• [-hI^{u'Y^) + e //^(uf ^) - Ill'^iuf^) - 1] ,(49) 
which in general differs from w^^^ = u^^^ + fw'-^''. From 



,(2) 



,(2) 4_ £-,,.(2) 



> it is possible in principle to solve 
for m(2) ^ from which (5^^) is similarly constructed, and we 
iterate henceforth via the prescription 

1 + (5("n*,^) _ 
1 + ^ 

- h{t) /^(ui")) + h'it) //^(uf ^) - iii''[u^Y^),{m 

that we use to integrate u;("+^). A consistent solution is 
then found in the limit 



£ = 7r-\ * 



n oo. 



(51) 



Note that after each iteration, the condition Vq xio^") = 
(arising from the Lagrange equations in the exact case) is 
preserved at all times. Finding a self-consistent solution 
to the dynamics by successive iterations is tantamount to 
minimization of the difference of the particle orbits 



AF{t,X) — Fsclf-consistont(^, ^) — F{t,X), 



(52) 



that we estimate at each iteration through the differential 
equation 



= (aa + 2d2) — AF(")+ad2 AjT^H (53) 



da 



da2 



This is however a non-trivial numerical problem that is 
beyond the scope of this article and we plan to pursue 
this possibility in a future work. 
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5. Comparison with other work 

We have calculated a unique and well-defined relation 
( p^ ) for cosmic density and velocity fields that is a self- 
consistent solution in the second-order model. Given the 
second-order ansatz for the particle orbits (||), the Lagrange- 
Newton system determines the functions of time /i (t) and 
/2(i) for generic initial conditions and it constrains the 
perturbation potentials S^^\ 5*^^^ in terms of the peculiar- 
velocity potential S* at to- Consequently, the orbits are de- 
termined at all times until shell-crossing, and the density 
is obtained from the continuity equation by evaluating the 
Jacobian of the transformation X ^ q. 

For simplicity, we have discussed the formulae on the 
assumption that the velocity and acceleration fields are 
parallel. In Sect. 3, it is shown that this assumption is a 
sufficient condition but by no means necessary, and in gen- 
eral it is sufficient to adopt the condition of parallelism of 
the gradients of the fields, which is a weaker requirement. 

The self-consistency of the solution for the density field 
follows from the fact that the orbits employed are second- 
order solutions in perturbation theory and not any given 
ansatz. Our result ( p^ ) generalizes previous results in the 
literature in the following respect. Nusser et al. (1991) 
obtained two relations for the density constrast in the ZA, 

6, {t, q) = -D I^'ivi) + 11^ {v,) - Ill^'iv,), (54) 



and 

5,{t,q)^-Dl''{vi), 

where a scaled peculiar velocity is used. 



(55) 



(56) 



and Sc is derived from the continuity equation and Se from 
Euler's equation. The difference between ( |5^ ) and ( p5| ) 
is second order and is therefore within the order of the 
approximation: 



Sc{t, q) - 6,{t, q) = 11^ {v,) - III^{vi). 



(57) 



Thus, in the ZA, only the first term on the RHS of ( p4| ) 
is significant and consistent with conservation of momen- 
tum. In the second-order result (22), we are able to retain 
two more terms and this relation has a similar form as 
(|5^), where D is replaced by h. This form of the over- 
density field satisfies automatically Euler's equation to 
second-order, since the orbits are obtained solving the 
Lagrange-Newton system to this order, and therefore the 
II{vi) and III{vi) terms are a legitimate part of the so- 
lution and not higher-order 'error' terms as in (^). 

One does not obtain self-consistent solutions in mod- 
els where an arbitrary perturbative parameter, such as 
D{t), is singled out, which results in a hybrid Eulerian- 
Lagrangian scheme. The chief problem inherent to such 
approaches is that the time dependence is set ab initio on 
the orbits and it is not determined by the equations of mo- 
tion as is done in the solutions we have used. This type of 



ansatz can be employed to good effect in a least-square fit 
for the dynamics of the particles as a two-boundary condi- 
tion problem (see for example Susperregi & Binney 1994), 
but in perturbation theory they are inconsistent and of 
little predictive power due to their lack of convergence as 
discussed in the Introduction. In such approaches, both 
the continuity equation and the equation of motion must 
be truncated to the given perturbation order and, conse- 
quently, the corresponding solutions 5c and 5^ in general 
differ. Therefore, from the analysis of the previous section, 
we can also conclude that 5c and 5c do not only differ by a 
higher order in the approximation, but failure to solve for 
the orbits perturbatively is conducive to discrepancies in 
the solutions at the same order. Given a parametrization 
of the orbits such as (||) , the second-order solution for the 
continuity density is given by (Gramann 1993b) 



5c = -D I^{Vr) + D'' 11^ {vi), 



(58) 



As shown by Gramann (1993b) the second order term for 
5c differs by a factor of 4/7 from that of (58). This incon- 
sistency at second order arises as a consequence of using 
second-order orbits that do not correspond to the actual 
second-order perturbation theory solutions. In contrast to 
the solution (^), (^ only accounts for quadratic cou- 
plings in Ui. Considering that at first order we already 
have the relation (|2^), it is apparent that (^) is not an 
accurate generalization of the Lagrangian linear theory. 

In the scheme discussed here, initial conditions are set 
aX t — tt), and hence we can try to establish a connection 
with the reconstruction model (^) by specifying these at 
t = 0. However, the normalization chosen for a{t) impedes 
to set ^0 identically equal to zero. A choice of an arbitrarily 
small to requires keeping only the growing solutions in the 
expression for h{t), entirely neglecting the presence of de- 
caying modes. We have to further restrict ourselves to the 
leading term (which is a particular solution of the second- 
order differential equation) and neglect the homogeneous 
second-order solution (which is growing). The leading or- 
der in h{t) is, at both first and second order (for model 
I), 



3 , 



1/2 



(59) 



and thus, using definition (|56| ) in (^ 

5^ -{had) /^(w,) + {haaYlI^{vi) 



j), we get 

- {haaf III^{vi) 



« -(a - 1) I^{v,) + {a- 1)2 //^(z;,) - (a - 1)^ III^{vi). 

(60) 

With the normalization chosen for a(t), ( pO| ) is in reason- 
able agreement with ( p8| ) at early times. At late times 
however, the presence of the III^{vi) term in ( |60| ) dimin- 
ishes the departures from linearity brought about by the 
quadratic term whereas this effect is not accounted for in 



It must be emphasized that ( p^ ) is strictly only valid 
so long as the inversion H{t, q) is well-defined. Indeed, 
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the expression given for 5{t, H{t, q)) provides the Eulerian 
values of the density contrast and is given as a function 
of Eulerian positions up to orbit-crossing time. At later 
times the inversion map becomes singular and this result 
does not hold. 

By integrating the Poisson equation in Eulerian coor- 
dinates using S{t, H{t, q) it is concluded (e.g. Gramann 
1993a, eq. (2.20)) that the presence of the quadratic term 
11^ (vi) in the density- velocity relation induces a depar- 
ture from parallelism between peculiar- velocity and -accel- 
eration vectors. The derivation of this relation under the 
assumption of parallelism proves this conclusion is not cor- 
rect. Certainly, one would be able to obtain agreement 
with solutions of Poisson's equation if a self-consistent 
density contrast were known. It is manifest however that 
this is far from being the case in these models. Differences 
in the solutions within the order of the approximation, as 
it is the case with Sc and Se in Gramann's analysis, imply 
that any conclusion drawn from Poisson's equation ap- 
plied to the reconstructed density is not justified as long 
as approximations are involved. 

We finally note that an expansion of the Jacobian to a 
given order (after the solutions arc obtained) mirrors for- 
mally the mathematical fact that, in principle, we cannot 
believe terms of higher order than the perturbative or- 
der of the trajectories. However, the powerful properties 
of a Lagrangian approach mainly rely on the existence 
of an exact integral for the density for any given trajec- 
tory field, which allows using these approximations until 
shell-crossing. At this epoch the contributions from sec- 
ond and third invariants is of the same order as that from 
the first invariant, a fact which underlies the success of the 
ZA. Whether this extrapolation until shell-crossing is jus- 
tified at all is another issue. The successes of Lagrangian 
perturbation methods as summarized by Bouchet et al. 
(1995) and the success of following trajectories until shell- 
crossing and even beyond (Melott et al. 1995) strongly 
argue in favour of this extrapolation. 

6. Conclusions 

In this paper we have presented a relation between cosmic 
density and velocity fields following a purely Lagrangian 
derivation (Eq. (p^) ) that is self-consistent in the second- 
order model. The advantage of working in the Lagrangian 
picture is that the over-density field obtained from the 
continuity equation is self-consistent and automatically 
satisfies conservation of momentum. For simplicity, one 
can adopt the condition ( |l9| ) of parallelism of the fields 
u and w, satisfied at first order and by a large class of 
irrotational solutions in second-order Lagrangian pertur- 
bation theory, and the result (^) holds accurately at both 
orders. However, as pointed out in Sect. 3, we note that 
( p^ ) is not a sine qua non condition but a weaker condition 
of parallelism of the gradients of the fields would yield the 



same result. It is shown that this condition is well satis- 
fied for generic initial conditions. In each case, the function 
h{t) is given explicitly in Appendix B.l and Appendix B.2 
respectively. Its range of validity is limited by the epoch 
when the inversion map H{t, q) becomes singular, i.e. up 
to orbit-crossing time, while its principal range of valid- 
ity should strictly be estimated by |pij(ic)| << a(^c), 
where p{t,X) — a{t) [F{t,X) — X] is the displacement 
field with respect to the Hubble- flow. This condition is 
probably very conservative in view of the success of these 
approximations if followed up to shell-crossing. 

By specifying the initial data at i « and keeping 
only the leading terms in the solutions, we are able to re- 
cover the leading terms of previous reconstruction models 
as shown in Sect. 5. Furthermore, the time dependence 
of the orbits in our model is established by solving the 
Lagrange-Newton system at each order. This leads to a 
unique and consistent result for the Eulerian density con- 
strast in terms of the peculiar-velocities. 

We finally wish to show some skepticism about the 
usefulness of the reconstruction formulae given. We want 
them to be applicable to present day observational data. 
However, we know confidently that shell-crossing has oc- 
curred, i.e. the interesting regime is no longer covered. 
Instead one would have to investigate a relationship be- 
tween smoothed cosmic over-density and peculiar-velocity 
fields, where the smoothing window scale has to ensure the 
absence of vorticity and multi-stream systems in the av- 
erage flow. Given the present work we should ask whether 
there is a formally proper way to incorporate smooth- 
ing into a consistent Lagrangian reconstruction formal- 
ism. With regard to this point, we think that solutions of 
the Lagrange-Newton system may possibly provide rea- 
sonable tools for reconstructing average fields, since aver- 
ages in Eulerian space of any vector function A(t, H{t, q)) 
are straightforward to calculate: {A)q oc J d?q A{t,X) — 
J (PX J{t,X) A{t,X). However, we still encounter the 
fundamental problem of averaging over multi-stream flows, 
which may not be properly handled in the framework of 
the Lagrange-Newton system. Rather, such a description 
has to be based on an approximation of the Vlasov-Poisson 
system. This major improvement of the methods lies well 
beyond the scope of present reconstruction techniques and 
requires the construction of Vlasov-Poisson type approxi- 
mations. 
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A.l. First- order model I 

The first-order particle orbits are 
Fit,X) = X + hit) to VxSiX), 



(Al) 



where the perturbation potential is related to the initial 
gravitational potential through the relation 

VxS = -tl Vx(i3- 



,i„(„l/2_^-l/2). 

Thus, we obtain: 



(A4) 

(A5) 
(A6) 



wto^ a^^t^^ VxS. 

From this we derive the well-known result for the Zel'dovich 
approximation used in the Sachs- Wolfe effect that 

V 



since 



const, 



w = 



at. 



2 V 



(A8) 







'ait) m' 






w{t,H) 




At) CW. 




yxs^^\x)_ 



where 
3 



a = to a 



f3 







,3/2 



,1/2 



C---^{l-a-y')-hl. 

Using ( A14 ), we invert (|^) in the form 

Hit,q) = q - (1 - a^'^/^)^^[miit) u + m2it) w] 

where 



3 7 5 
(A2) mi ^ -to i-a^^^ - -a- 



1/2 



7 

20^ 



-3/2 



9 

(A18) 
(A19) 
(A20) 



The inversion of (Al) is then: 

Hit,q) ^ q-hit) uit,H) = q-a~^'^ to hit) i?),(A3) 
and therefore hit) is given by: 
3 



Therefore, given the condition of parallelism (|19|), we have 
that ( |A19D reduces to (|o|) with 

(A22) 



hit) = (1 - a-^'^)-\miit)+t-^m2it)] 



= ^ to (1 - a-^/2)-i(ai/2 _ „-i/2 _ „^3 ^ ^-4)^ (^23) 



A. 3. Second- order model II 



Following the same procedure as in Appendix A. 2, we 

(A7) 

write 



" uit,H) ' 




'ait) Pit)' 




'S^xS^'HX)' 


wit,H) 








.Va-5(2)(X)_ 



where 



A. 2. Second-order model I ^ 
Eqs. (|l7|)([l8|) in terms of the perturbation potentials read p 

(A9) 



2,^-l(„l/2+3^-2)^ 



9 

25 



2 3/2_|_^„l/2_^l^-l^ 



Lt-i(_±a3/2 



(A24) 



(A25) 

1^-7/2) (A26) 
70 4 ''^ ' 



V = 



UoHa-'-a-y% 



^xS^'HX) 
S/xS^^Hx) 



7 5 35 

V = to^ 

^ = ^° (-14+10" +35" 
We write the perturbation potentials in terms of u and w a = 28 to a 
by inverting (|A9|): 



(AlO) 
(All) 
(A12) 

(A13) where 



3 2/ 10 29 1 43 7/2 5 
^=25*0 (--+20« +35«^ ' - r ) 



Therefore, inversion of (A24) leads to 

= 7W 



a'it) P'it) 

v'{t) at) 



uit,H) 
wit,H) 



(A27) 
(A28) 

(A29) 
(A30) 



/3' 



'VxS^'Hx)' 




■ a'it) P'it) ' 




' uit,H) ' 


.Vx^(2)(^). 




.i{t) Cit)_ 




wit, H) 



(A14) 



where 

a' = (l-a-V2)-4^ (5^-1/2 



21 
10' 



-3/2 



P'^il-a-yT'tli-la+^^-la'^/% 
V'^lil-a-y^'toa-'/^ 



(A15) 
(A16) 7 
(A17) 



y 2 / 

"5*0^7" '60 
7^' ^2toi-a-' + a-y') 
C^tl (2ai/2 + 3a-2), 



2^3/2^29^1/2^1 _i_43^_2^1 -7/2 



70 



(A31) 

(A32) 
(A33) 



/ 24 21 2 45 , 9 n/2 27 7\-\,„,, 

a^'^ a^^H a^^ a^^'^ a^^] .(A34) 

V 175 25 28 25 100 / ^ ^ 



10 



M. Susperregi & T. Buchert: Cosmic density and velocity fields in Lagrangian perturbation theory 



Appendix B: Jacobian 
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The Jacobian det{d Hi /dqj) is given by: 

J^^l-g, lE{u,) - 52 I^'iw,) + gl II^{u{) + gl II{w{) 
+31 32 ^"'{u^ I w,) - gl III^{u,) - gl III^{w,) 
-gl 32 ///'^(m^ \w,)~ 51 gl in"^{u, \ Wj), 



where 



II' {Ui \Wj) = -{Ui,iWj^j - Ui^jWj,i), 

III^{ui) = det(itjj) 



(Bl) 

(B2) 
(B3) 

(B4) 
(B5) 



IIl"^{u^ \Wj) = 
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+"3,2^1, 1^2, 3 - W3,lU'l, 2^2,3 " "2,3^1,2^3,1 + "2, 21^1,3^3,1 
+Ml,3W2,2t«3,l - Ul, 2^2,3^3,1 + W2,3Wl,lU'3,2 " ^2, lU^l, 3^3,2 , 
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(B7) 

and the functions 51 (t) and 52 (t) are model-dependent. 
For instance, for the case of model I of Sect. 3, 



5i(t) = (l-a-7/2)-imi(i), 
52(<) = (l-a-7/2)-i 



(B8) 
(B9) 



where mi(t) and m2{t) are defined in Appendix A. 2. The 
assumption of paralleUsm (19) adopted, is simplified 
by making the following substitutions in (p3l|): 



w ^0 

.91 W ^ Ht). 



(BIO) 



where 

h[t)^g^{t) + r^g2{t). 



(Bll) 



Hence, all terms containing Wi are set equal to zero, 
and thus (p31|) simply reads 



J'^ = \-h I^{ui) + II^{u,) - Ill^iui). (B12) 



